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0- . Abstract 

o ■ 

I/"") ■ The purpose of this work is to compare the dynamics of arrays of Josephson junc- 

tions in presence of magnetic field in two different frameworks: the so called XY 
frustrated model with no self inductance and an approach that takes into account the 
screening currents (considering self inductances only) . We show that while for a range 
of parameters the simpler model is sufficiently accurate, in a region of the parameter 

; 1 Introduction 

Arrays of Josephson junctions have been proposed more than two decades ago to enhance 
the emission of microwave [|]] . In fact it is well known that the power available from a single 
junctions is not enough for many practical applications and therefore the achievement of 
coherent motion of arrays of junctions is an important issue for device applications [§]. Apart 
from the applications, two-dimensional arrays have been also investigated as an interesting 
nonlinear system both experimentally and theoretically [|3|, |4|. In the study of such arrays 
two classes of models have been proposed to describe the dynamics: 

1) Models that neglect the effect of the screening currents. Those model are often called 
"Uniformly Frustrated XY" j5], ||, [7[] because the Hamiltonian for this model is similar 
to that of a square flat lattice of spins in magnetic field (for that reason it is also called 
Spin Glass model). 

2) Models that include, to some extent, the effect of the screening currents. The simplest 
version of this model assumes that the effect of the screening is to shield the magnetic 
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field only in the two cells adjacent to each branch where a current is flowing (in the 
following this model will be termed the NS model after Nakajima and Sawada that 
have introduced it |8[]). More generally one should consider the effect of the screening 
in all cells (full mutual inductance approach |10||). 

To neglect screening currents greatly simplifies the problem. Moreover, the formal simi- 
larity with other well known systems allows to take advantage of an accumulated experience 
in those contexts. Model (2) has been used less extensively also because, in its complete form 
(that includes all mutual inductances), requires complex numerical routines and the result- 
ing equations are difficult to handle analytically []9], |1(|] . In spite of the simplicity of model 
(1), it has been able to explain most of the experimental observations on two-dimensional 
arrays, for example giant Shapiro steps ||llfl . It is therefore of interest to establish the limits 
of validity of this approach, that should be taken into account when interpreting the exper- 
imental data. Moreover, to develop coherent sources based on large two-dimensional arrays 
of Josephson junctions, better performances can be obtained if the junctions in the array 
have larger critical currents; in fact the emitted power of a single junction is proportional 
to the square of the critical current, while the emitted power of an ideal coherent array will 
oc N 2 Iq, where N is the number of locked junctions in the array Hi 2(1 . This, as will be shown 
in the following, prevents for practical purposes to use arrays whose coupling is so strong to 
avoid large screening currents and therefore the model (2) seems most suitable to describe 
arrays usable as microwave sources. 

The purpose of this work is to explore the parameter space looking for the limits of 
validity of model (1). We will show that, for certain ranges of the parameters and at least 
for the simplest arrays, the solutions of model (1) are really different from the solutions of 
model (2). To fully understand the differences between the two approaches we will present a 
step by step derivation of the dynamical equations. As a tutorial example we will start with 
a discussion of a simple dc SQUID (Superconducting QUantum Interferometer) containing 
only two junctions. This well known case will lead us to some general considerations on the 
role of the screening currents. 

It will be shown that the XY equations cannot be easily derived from NS equations 
taking some limits of the parameters. It will be further shown numerically that there is a 
region of the parameter space where it is not possible to assume that the screening current 
is negligible. Next, we will consider the most elementary cell of a two-dimensional array, i.e. 
a square cell with junctions (and inductances) on each branch |]13(] . We will again derive 
the equations and we will numerically show that the screening current are not negligible for 
some ranges of the parameters. A comparison with the actual two-dimensional arrays will 
be performed for some cases to check that the results for the single cell are reliable. Finally, 
the macroscopic dynamics in the two regimes (with negligible and not negligible screening 
currents) will also be shown. 



2 The Squid model 



A dc SQUID is a superconducting loop interrupted by two junctions [|T2|| . Since there is a 



single superconducting loop the flux quantization, neglecting screening currents, is given by: 
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0! - 2 = 2nf (1) 

where 0i and 2 are the gauge invariant phases across the two junctions, / = $ a /$o is 
the frustration, or the ratio between the applied flux $ a and the quantum flux $ = h/2e. 
Assuming the RSJ model for the junctions [|12] the equation of motion of the phase difference 



is given by the current balance across a branch. In normalized units it reads: 

01 + ct0i + sin 0i = 7, (2) 

02 + a02 + sin0 2 = 7. (3) 

Here time is normalized with respect to the inverse of the Josephson frequency u>j = 
27rJ /C$o (C is the capacitance of the junction and I is the maximum Josephson cur- 
rent), a = ($o/2vrCi? 2 / ) 1//2 is the damping term (R is the shunt resistance), 7 is the bias 
current normalized to Jo- Here and in the following we will assume the parameters of the 
junctions to be identical. 

Taking the sum of Eq.s (2,3) and inserting Eq. (1) the equation of motion for 0i is readily 
obtained: 

7 = 0i + «0i + - [sin 0i + sin(0i - 2irf)] . (4) 

Under the approximation of negligible screening current the SQUID is described by Eq. 
(4). Although it is possible to recognize some characteristics of the SQUID, in some respect 
Eq. (4) is qualitatively incomplete. To better recognize the difference it is useful to follow 
the complete derivation of the equation considering also the screening current [|14j| , i.e., 
modifying Eq. (1) as follows: 

2ttL 

0i - 2 = 2nf + —I, (5) 
<±>o 

(here L is the inductance of the superconducting loop and I s is the screening current in the 
loop). The current balance now is modified by the presence of the screening current: 

01 + a0i + sin0i = 7-7^ (6) 

02 + a02 + sin0 2 = 7 + t^- (7) 

Inserting Eq. (5) will lead to two coupled equations rather than one is the SQUID 
parameter = 2tiLI /& ) 

... 1 2tt 

01 + «0i + sin 0i = 7 - —(02 - 0i) + — /, (8) 

Pi Pi 
... 1 2tt 

02 + «02 + sin 02 = 7 + -5- (02 - 0i) — w~ f ■ (9) 



Pi TX ' A 



The most obvious difference between the Eq.s (8,9) and the XY model is the increased 
number of equations: while in the XY approach the dynamics of the junctions are essentially 
identical [apart from an additive constant, see Eq. (1)] in the NS context the two dynamics 
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are governed by two coupled differential equations. This remarkable difference is a conse- 
quence of the fact that the screening current is a dynamical variable itself, and the fluxoid 
quantization can be used to eliminate it from the equations in the NS approach, whereas 
the quantization rule can be used to eliminate a phase variable in the XY approach. A 
common characteristic of the two sets of equations is that for / = a solution is the single 
junction free running solution. The difference is that while for the XY system this is the only 
solution, the NS system can allow for other types of solutions, at least for some parameters 
values. The NS model allows the presence of "propagating solutions" , in the sense that an 
excitation in a junction can propagate to the other, giving rise to the so called "beating 
solutions" | ]15| , |16| j . As a consequence the NS approach will reveal a much richer dynamics 



than the XY approach; for instance, it can show hysteresis also in the limit of negligible 
capacitance (in this limit the second derivative term disappears), while the XY equation 
cannot. Another remarkable difference is that for / = 1/2 the XY equation becomes linear, 
while the NS set of equations still retains its nonlinear terms. 

To get a deeper insight into the difference between the two approaches we have numer- 
ically evaluated I s integrating the Eqs. (8,9). To estimate the importance of I s in the 
dynamics we have plotted the maximum of its absolute value: 

J™ = max\I a /I \. (10) 

Initial conditions are always chosen as 0i = 02 = sin -1 7 if 7 < 1 and zero otherwise, 
moreover fa was set to zero. We have used two numerical methods (a simple 4'th order 
Runge Kutta and a more elaborate Bulirsch Stoer) obtaining consistent results. In Fig. 1 
we show I™ as a function of the bias current 7. The two curves are obtained slowly increasing 
(diamonds) or decreasing (crosses) the bias current. For both curves there is a region where 
the screening current can reach significant values. To check that I™ is really a significant 
test to discriminate where solutions of the NS equations are different from those of the XY 
model, we have plotted the microscopic dynamics in Fig. 2 for three points, two for zones 
of the bias in which I™ is not negligible and one for higher bias value. It is evident that in 
the first two cases the difference 02 — 0i is not just a constant as assumed by the XY model 
[see Eq. (1)], while this is indeed a good approximation when 7™ is negligible. 



3 The cell model 

As a first step toward the study of two-dimensional arrays we will consider in this section 
the most elementary cell of a square array (see Fig. 3). This cell has already been used to 
infer the properties of two-dimensional arrays (with rf bias) by Sohn and Octavio ||13 |, a 
row of such cells has been considered to investigate the stability of Josephson array ||17|1 . In 
the XY approach the equation for the flux quantization reads (see Fig. 3 for notation): 

V x - H 2 - V 2 + H x = 2nf (11) 



Moreover, for the symmetry of the system also the following identity holds []T^, |17 



Ht = -H 2 = H. (12) 
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The current balance in the four nodes can be written as (Jx denotes the current flowing 
through the RSJ elements of the junction X): 



7 = J H - Jvi (13) 

7 = ~Jn-Jv % (14) 

7 = ~Jh ~ Jv 2 (15) 

7 = J H - Jvi (16) 

Taking the sum and difference of the independent equations the dynamics of the cell in 



the XY model is governed by the following two equations: [13 



S + aS + 2an^-S^cos(-D\ = -2 7 (17) 
D + aD + sin ( ~d) cos ( -s) + sin ( -D -irf) = (18) 



where S = V x + V 2 and D = V x - V 2 . 

For the NS approach, as usual, we have to consider also the screening current, while for 
the same symmetry reasons Eq. (12) still holds. In conclusion Eq. (11) becomes: 

V 1 -V 2 + 2H = 2itf - (19) 

and the equations of motion of the cell are: 

•■ 1 2tt 

V 1 + aV 1 + smV 1 = -7 + - {V x - V 2 + 2H) - — f (20) 

Pi Pi 

V 2 + aV 2 + smV 2 = -7 - - (Vt - V 2 + 2H) + ^/ (21) 

pi pi 

H + aH + sinH = \ {V x - V 2 - 2H) - ^/ (22) 

Pi Pi 

For comparison we write also Eq.s (20-22) in terms of the variables S and D; in this case 
the Eq.s yield: 

S + aS + 2sm(^S^ cos (~Dj = -2 7 (23) 

D + aD-2cosQs)sinQ£>) = j(D + 2H)-j-f (24) 

H + aH + sinH = l(D + 2H)-^f (25) 

Pi Pi 

Eq. (23) and Eq. (17) are identical, but Eq. (18) is not easily recognized as an ap- 
proximations of Eq.s (24,25). Therefore the considerations done for the SQUID apply also 
here: the two approaches differ in the number of equations and it seems difficult to predict a 
priori for which parameter values the screening current is negligible. Rather, we numerically 
integrate Eqs. (20-22) and in Fig. 4 we show the dependence of I™ on the parameter /3i as 
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a function of the bias current for two different values of the damping coefficient a; in Fig. 5 
is shown the behavior of the screening current as a function of the frustration and the bias 
current, for the same values of the damping coefficient. The results shown in these figures 
can be summarized as follows: 

1) For 7 >> 1 the screening current decreases. This corresponds to the observation that 
(in the overdamped limit) for high bias current the solutions are nearly sinusoidal (plus 
a constant slope) and with a constant phase-shift in presence of magnetic field []17[J . 
This state corresponds to negligible screening current. 

2) For high inductance values (for instance when fii > 2 for a = 0.25 and / = 0.5) the 
screening current is also negligible. On the contrary, for the same values of a and / and 
for Pi values down to 0.1 (a fairly low value for practical arrays [[TJ|]), the screening 
current increases, up to values where it is not negligible; however for Pi < 0.1 this 
region shrinks. 

3) The capacitance and the dissipation play an important role in controlling the screening 
current effect. Significant screening current are observed for a < 1, i.e., for relatively 
high values of the capacitance or for low values of the dissipation. 

It is also important to notice that, as expected, hysteresis plays an important role. In 
the figures shown so far we have always used as initial conditions all the phases equal to 
zero and at rest. As already shown in Fig. 1 a different choice of initial conditions can lead 
to a change of the results; This was observed, rather than changing the initial conditions on 
the phases, simply sweeping the current bias starting from values different from zero. In the 
cell case simulations shown that there are solutions in the underdamped limit that exhibit 
a significant contribution to the screening current also in absence of magnetic field. 

As we have pointed out in the Introduction, two-dimensional arrays are of interest if a 
large number of junctions can be locked together. Even if we believe that the simple cell 
studied in the previous section offers the distinct advantage of simplicity and illustrates the 
basic mechanism, it is quite natural to ask if it can also furnish quantitative predictions on 
larger arrays. To check that the results obtained so far are not crucially dependent on the 
fact that we are considering only an elementary cell, we have investigated the behavior of 
larger arrays in few cases, using the NS equations [|8|]. In Fig. 6 it is shown the behavior of 
a 10 x 10 array where zones of screening current comparable with the elementary cell can be 
clearly seen; in that zones the overall dynamics can be directly compared with that of the 
elementary cell (see also discussion below) . These zones corresponds to a large range of the 
parameters, however, besides these zones in the low damping case (a = 0.25, see Fig. 6a), 
a region of very large maximum screening current is shown: this region corresponds to the 
penetration of static fluxons in the array and does not have any correspondence in the cell 
case; further studies of this region will be carried out in the future. 



6 



4 Microscopic dynamics 



So far we have dealt only with the "macroscopic" quantity I™. Although we claim that 
this is sufficient to discriminate between states were model (1) and model (2) do not give 
consistent results (because the set of equations are the same if I s is negligible) it is interesting 
to investigate in which sense the dynamics differ when the screening current is not negligible. 
For sake of simplicity we will concentrate on the microscopic dynamics in the elementary 
square cell. 

We have systematically found that when the screening current is not negligible there 
appear solutions that are close to the well known "beating solutions" for SQUIDs ||15 , |16 



i.e., to solutions that correspond to the entry and propagation of ffuxons across the cell, 
or to the motion of fluxons across the array. To recognize this we have shown in Fig. 7 
the maximum screening current as a function of the voltage behavior together with the IV 
characteristic. It is evident that a large screening current is associated to a resonant step that 
is reminiscent of the Fiske steps in long Josephson junctions. In fact, in the long Josephson 
junction language resonant structures are named Fiske step if they occur in presence of 



magnetic field, and zero-field steps if they occur in absence of magnetic field [12 



It is also evident that there is a sharp change in the behavior of I™ when the systems 
switches from the resonant step to another solution. To better understand the nature of 
this transition we have plotted in Fig. 8 the voltage across the two vertical junctions for 
two values of the bias current: just before the switch and just after. The difference between 
the two dynamics is that while on the resonant step there is a large voltage pulse, after the 
switch the dynamics is much more uniform with a smaller modulation of the voltage. We 
claim that the transition is similar to that observed for one- dimensional arrays or for long 
Josephson junctions ||19|, |2~0 . 



A comparison between Fig.s 1,2 and Fig.s 7,8 strongly indicates that the dynamics are 
very similar to that observed for the simple SQUID: in fact as in the SQUID case resonant 
states appear in general together with high values of screening currents and similarly the 
dynamics appears to be uniform where screening currents are negligible. 

The different nature of the two solutions can be clearly seen in Fig. 9, where we show 
a plot of the time delay across the two vertical junctions of two peaks of voltage. Indeed, 
while on the resonant step this time delay roughly tends to a constant (corresponding to 
the maximum speed of propagation of a signal across the system), in the state where the 
screening currents are negligible the speed increases with current. This suggests that in the 
latter mode there is no actual signal propagating across the cell, but rather a modulated 
solution similar to that obtained in Ref. ||17|| for overdamped junctions in the high bias 
limit. It should be noticed that the very concept of a propagating fluxon is not applicable to 
the solution arising from the XY model: If we adopt as a working definition of propagating 
fluxon an excitation of roughly one flux quantum moving from a cell to another, then to 
neglect screening currents prevents the appearance of such solutions because the magnetic 
field is supposed to be uniform across the array. In this language the XY model does not 
allow signal propagation, so we conclude that the XY model can be an approximation for 
the NS equations (as well as any model including screening current effects through mutual 
inductances) only if the latter do not carry a propagating solution. 
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Since the uniform regime occurs for higher bias, an interesting question is if the border 
between those two regions corresponds to the border between the regions of flux motion 
and uniform solutions devised by Lachenmann et al. with the LTSEM (Low Temperature 
Scanning Electron Microscopy) technique ||21|| . When comparing the thresholds obtained 
numerically with those observed experimentally it is quite arbitrary to decide how small 
the screening current should be to lead to the disappearance of the signal measured by the 
LTSEM, therefore a detailed quantitative comparison is not possible on the basis of the 
results shown here. However, the measured threshold for an array 10 x 10 critically damped 
(a ~ 1) is 7 ~ 4 |]21[|, a region where our simulations predict that the screening current is 



small: I™ ~ 0.16 or roughly 4% of the bias current at most, i.e., for / = 0.5. 

It is nevertheless possible to notice that both the screening current and the LTSEM 



signal indicates an homogeneous solution for high bias current H21[] . If the hypothesis that 



the two thresholds are related is true, it is possible to go further and to speculate that the 
coherent emission from Josephson junctions arrays is related to the presence of non negligible 
screening currents. 



5 Conclusions 

We have proven that the choice of the more appropriated model for Josephson junction 
arrays depends on the parameters, especially on the bias current and the SQUID parameter 
It is possible in fact to show numerically that for certain regions of the parameters 
the contribution to flux quantization arising from the screening current, neglected in the 
so-called XY model, is important, and this can happen also for relatively small values of the 
inductance of the loop. We have traced the origin of this to the structure of the equations, 
noticing that the screening current is a dynamical variable itself whose value cannot be easily 
predicted a 'priori, and we were able to investigate its behavior only numerically. We have 
also found that the presence of a non negligible screening current prevents the occurrence 
of uniform solutions to set in and rather induces solutions that are known, in the context 
of SQUID's, as "beating solutions" pT5 , Op . Finally, there are reasons to believe that 



there might be a connection between the presence of significant screening currents and the 



phenomenon of "row switching" observed with the technique of LTSEM [|21[ . 

It is worth to recall that we have a) neglected any spread of the parameters, b) retained 
self inductances only and neglected mutual inductances. 
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Figure Captions 



Maximum screening current versus bias current for a SQUID. Parameters of the sim- 
ulations are: a = 0.25, Pi = 1.0, / = 0.5. 

Dynamics of the phase difference of the two junctions of the SQUID for 7 = 0.18 (solid 
line), 7 = 0.70 (short dotted line), and 7 = 1.0 (dotted line). Other parameters are 
the same as in Fig. 1. 

Schematic of the elementary cell of the two-dimensional square array. Crosses denote 
the Josephson elements. 

3D plot of the maximum of the screening current I™ versus the bias current 7 and the 
SQUID parameter Pi for the elementary cell for (a) a = 0.25 and (b) a — 1.0. The 
frustration is set to / = 0.5. 

3D plot of the maximum of the screening current I™ versus the bias current 7 and the 
frustration / for the elementary cell for (a) a = 0.25 and (b) a = 1.0. The SQUID 
parameter is set to Pi = 1.0. 

Maximum screening current as a function of 7 and / for an array 10 x 10 with Pi — 1. 
Parameters of the simulations are: (a) a = 0.25, (b) a = 1.0. 

Voltage (dotted line) and of I™ (solid line) as a function of the bias current in the 
square cell. Parameters of the simulations are: a = 0.25, Pi = 1, / = 0.5. 

Dynamics of the voltage across the vertical junctions in the square cell for 7 = 0.84 
(solid line and dotted line) and 7 = 0.86 (dotted lines with shorter length). Parameters 
of the simulations are the same as Fig. 7. 

Time delay between two peaks of the voltage in two adjacent vertical junctions for 
Pi = 1 (+) and Pi = 0.5 (stars). Parameters of the simulations are: a = 0.25, / = 0.5. 
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